The Journal of Chemical Physics

Biomolecular condensates, membrane-less organelles that arise from liquid-liquid phase separation (LLPS) of proteins and nucleic acids, play vital roles in cellular organization and regulation. Computational modeling is crucial for uncovering the molecular mechanisms behind LLPS; however, the fluctuations across a wide range of spatial scales and the inherently non-equilibrium nature of these systems make capturing their long-timescale dynamics particularly challenging. Here, we present a continuum dynamical density functional theory (DDFT) framework that captures the non-equilibrium dynamics of LLPS by integrating a physics-based statistical mechanical theory with key experimentally-derived parameters. Our model couples DDFT with a continuum free energy functional, incorporating two-body correlations between monomers and surface tension effects to determine binodal densities under phase coexistence. By solving the DDFT equations, we describe the time evolution of phase-separated domains, capturing key long-timescale processes such as droplet maturation, coalescence, and interface relaxation, phenomena that are difficult to probe using atomistic or mesoscale coarse-grained simulations. This implementation integrates experimental phase equilibrium data with molecular-scale descriptors, such as amino acid properties, to construct a quantitative link between molecular interactions and macroscopic phase behavior. However, the approach is generalizable, providing a foundation for self-contained molecular-to-continuum modeling bridge platforms.

Electron beam irradiation can cause damage to biological and organic samples, as determined via transmission electron microscopy (TEM). Cryo-electron microscopy (cryo-EM) significantly reduces such damage by quickly freezing the environmental water around organic molecules. However, there are multiple hypotheses about the mechanism of cryo-protection in cryo-EM. A lower temperature can cause less molecular dissociation in the first stage, or frozen water can have a "cage" effect by preventing the dissociated fragments from flying away. In this work, we used real-time time-dependent density functional theory (rt-TDDFT-MD) molecular dynamic simulations to study the related dynamics. We used our newly developed natural orbital branching (NOB) algorithm to describe the molecular dissociation process after the molecule is ionized. We found that despite the difference in surrounding water molecules at different temperatures, the initial dissociation process is similar. On the other hand, the dissociated fragments will fly away at room temperature, while they will remain in the same cage when frozen water is used. Our results provide direct support for the cage effect mechanism.

Building on a complex between a tubulin protofilament (PF) and a fragment of the Tau protein containing residues 169 to 367, we investigate the dynamics of the disordered elements of the system, namely the tubulin C-terminal tails (CTTs) and the Tau protein, using classical all-atom molecular dynamics simulations. Our results show that CTTs adopt a hook-like dynamic pattern on the bare PF while remaining highly mobile. The binding of Tau on the PF surface alters the dynamics of the I-CTTs in a sequence-dependent manner. While the repeat domains of Tau are mostly maintained on the PF by weak and strong binding patches with the tubulin cores, the Proline-Rich Region (PRR) relies on the wrapping phenomenon of I-CTTs to fuzzily stabilize its interaction with the PF. Our study thus provides a deep dive into the dynamic interplay between the Tau protein and the CTTs of microtubules, the latter being characterized extensively using a variety of disorder-adapted metrics. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=111 SRC="FIGDIR/small/721901v1_ufig1.gif" ALT="Figure 1"> View larger version (25K): org.highwire.dtl.DTLVardef@b3f985org.highwire.dtl.DTLVardef@1c2bf70org.highwire.dtl.DTLVardef@a66b95org.highwire.dtl.DTLVardef@1e138e0_HPS_FORMAT_FIGEXP M_FIG C_FIG

Biomolecular condensates are important contributors to the internal organization of the cell material. While initially described as liquid-like droplets, the term biomolecular condensates is now used to describe a diversity of condensed phase assemblies with material properties extending from low to high viscous liquids, gels, and even glasses. Because the material properties of condensates are determined by the intrinsic behaviour of their molecules, characterising such properties is integral to rationalising the molecular mechanisms that dictate their functions and roles in health and disease. Here, we apply and compare three distinct computational methods to measure the viscoelasticity of biomolecular condensates in molecular simulations. These methods are the shear stress relaxation modulus integration (SSRMI), the oscillatory shear (OS) technique, and the bead tracking (BT) method. We find that, although all of these methods provide consistent results for the viscosity of the condensates, the SSRMI and OS techniques outperform the BT method in terms of computational efficiency and statistical uncertainty. We, thus, apply the SSRMI and OS techniques for a set of 12 different protein/RNA systems using a sequence-dependent high-resolution coarse-grained model. Our results reveal a strong correlation between condensate viscosity and density, as well as with protein/RNA length and the number of stickers vs. spacers in the amino-acid protein sequence. Moreover, we couple the SSRMI and the OS technique to nonequilibrium molecular dynamics simulations that mimic the progressive liquid-to-gel transition of protein condensates due to the accumulation of inter-protein {beta}-sheets. We compare the behaviour of three different protein condensates--i.e., those formed by either hnRNPA1, FUS, or TDP-43 proteins--whose liquid-to-gel transitions are associated with the onset of amyotrophic lateral sclerosis and frontotemporal dementia. We find that both SSRMI and OS techniques successfully predict the transition from functional liquid-like behaviour to kinetically arrested states once the network of inter-protein {beta}-sheets has percolated through the condensates. Overall, our work provides a comparison of different modelling rheological techniques to assess the viscosity of biomolecular condensates, a critical magnitude that provides information on the behaviour of biomolecules inside condensates.

The DNA-binding protein from starved cells (Dps) compacts bacterial DNA into stress-protective condensates, yet the physical mechanisms underlying this process and the material properties of the resulting condensates remain poorly understood. Here, we combine coarse-grained Brownian dynamics simulations with Flory-Huggins polymer theory to elucidate the structural, dynamic, and thermodynamic principles governing Dps:DNA organization and condensate formation. The simulations, in which DNA is represented as bead-spring polymers and Dps as spherical particles, reveal that weak Dps:DNA attraction and low Dps concentrations produce extended, network-like morphologies, whereas stronger interactions and higher Dps concentrations drive compaction into dense globular condensates with suppressed DNA mobility and sub-diffusive dynamics. Complementary Flory-Huggins analysis identifies the corresponding thermodynamic regimes and shows how Dps:DNA affinity, DNA:DNA and Dps:Dps repulsion, and solvent quality determine the boundaries between homogeneous and phase-separated states. Together, the two approaches provide a unified microscopic and thermodynamic picture of condensate formation, bridging molecular interactions with emergent mesoscale structures. These results advance understanding of protein-nucleic-acid phase behavior and illustrate general principles governing biomolecular condensation in soft matter systems.

The spatial architecture and mechanical rigidity of polysomes are crucial determinants of translational efficiency and mRNA stability. In this study, we investigate the conformational statistics of an mRNA backbone decorated with high-density ribosomes at varying densities using large-scale, extensive molecular dynamics simulations based on the Kremer-Grest bead-spring model. To address the extreme spatial asymmetry between mRNA monomers and ribosomes, we used an efficient tree-based neighbour list algorithm, enabling the analysis of mRNA chains up to N = 4, 969. Our results demonstrate that the excluded volume of massive ribosomes induces a significant and robust expansion of the scaling exponent v from 0.59 to approximately 0.7. In the conformation of mRNA, this shift translates to a self-induced dimensional reduction from a three-dimensional random coil toward a stretched, a quasi-two-dimensional architecture at biologically relevant scales. Such a transition is further evidenced by a periodic "regain" of the bond-bond correlation function C(n) at ribosome attachment sites, indicating a geometric alignment absent in standard homopolymers. These findings reveal that the geometric crowding of ribosomes itself provides a robust physical prerequisite for the formation of higher-order polysome architectures, bridging the gap between polymer physics and structural properties of mRNA during translation.

Nuclear pores are protein assemblies inserted in the nuclear envelope of eukaryotic cells, acting as main gates for communication between nucleus and cytoplasm. So far, nuclear pores have been extensively studied to determine their structure and composition, yet their spatial organization and geometric arrangement on the nuclear surface are still poorly understood. Here, we analyze super-resolution images of the surface of Xenopus laevis oocyte nuclei during development, and characterize the arrangement of nuclear pores using tools commonly employed to study the atomic structural and topological features of soft matter. To interpret the experimental results, we hypothesize an effective interaction among nuclear pores and implemented it in extensive numerical simulations of octagonal clusters mimicking typical pore shapes. Thanks to our simple model, we find simulated spatial distributions of nuclear pores that are in excellent agreement with experiments, suggesting that an effective interaction among nuclear pores exists and could explain their geometrical arrangement. Furthermore, our results show that the statistical features of the geometric arrangement of nuclear pores do not depend on the type of pore-pore interaction, attractive or repulsive, but are mainly determined by the octagonal symmetry of each single pore. These results pave the way to further studies needed to determine the biological nature of pore-pore interactions.

Intermolecular spin relaxation by translational motion of spin pairs have been widely used to study properties of the biomolecules in liquids. Notably, solvent paramagnetic relaxation enhancement (sPRE) arising from paramagnetic cosolutes has gained attentions for various applications, including the structural refinement of intrinsically disordered proteins, cosolute-induced protein denaturation, and the characterization of residue-specific effective near-surface electrostatic potentials (ENS). Among these applications, the transverse sPRE rate known as {Gamma} 2 has been predominantly been interpreted empirically as being proportional to <r-6>norm. In this study, we present a rigorous theoretical interpretation of {Gamma} 2 that it is instead proportional to <r-4>norm and provide explicit formula for calculating <r-4>norm without any adjustable parameters. This interpretation is independent of the type or strength of interactions and can be broadly applied, including to the precise interpretation of ENS.

The antagonistic interaction between small RNAs (sRNAs) and messenger RNAs (mRNAs) constitutes a fundamental regulatory mechanism of gene expression in both prokaryotic and eukaryotic cells. However, the stochastic nature of transcription renders mean-field approximations inadequate for quantitative analysis of such systems. In the regime of strong sRNA-mRNA antagonism, we generalize the conventional probability-generating-function (PGF) framework and derive a novel approximate solution in the form of a generalized PGF, which can be analytically transformed into the time-dependent joint distribution of sRNA and mRNA via Laurent series expansion. The proposed approximation accurately captures the full stochastic dynamics across diverse systems exhibiting strong antagonism, while incorporating key biological features such as transcriptional burstiness, translation and sRNA recycling over the entire temporal range. Building on this analytical foundation, we further develop a generalized-PGF-based parameter-inference method that enables efficient and precise estimation of kinetic parameters, achieving inference speeds up to three orders of magnitude faster than traditional maximum-likelihood estimation approaches.

The ability to predict and understand the complex molecular motions occurring over diverse timescales ranging from picoseconds to seconds and even hours occurring in biological systems remains one of the largest challenges to chemical theory. Markov State Models (MSMs), which provide a memoryless description of the transitions between different states of a biochemical system, have provided numerous important physically transparent insights into biological function. However, constructing these models often necessitates performing extremely long molecular simulations to converge the rates. Here we show that by incorporating memory via the time-convolutionless generalized master equation (TCL-GME) one can build a theoretically transparent and physically intuitive memory-enriched model of biochemical processes with up to a three orders of magnitude reduction in the simulation data required while also providing a higher temporal resolution. We derive the conditions under which the TCL-GME provides a more efficient means to capture slow dynamics than MSMs and rigorously prove when the two provide equally valid and efficient descriptions of the slow configurational dynamics. We further introduce a simple averaging procedure that enables our TCL-GME approach to quickly converge and accurately predict long-time dynamics even when parameterized with noisy reference data arising from short trajectories. We illustrate the advantages of the TCL-GME using alanine dipeptide, the human argonaute complex, and FiP35 WW domain.

The phase separation of lipid bilayers, composed of mixtures of saturated and unsaturated lipids and cholesterol, is a topic of fundamental importance in membrane biophysics and cell biology. The formation of lipid domains, including liquid-disordered domains enriched in unsaturated lipids and liquid-ordered domains enriched in saturated lipids and cholesterol is believed to be essential to the function of many membrane proteins. Experiment, theory, and simulation have been used to develop a general understanding of the thermodynamic driving forces underlying phase separation in ternary and quaternary lipid mixtures. However, the kinetics of early events in lipid phase separation in the presence of transmembrane proteins remain relatively understudied. Using large-scale all-atom and coarse-grained simulations, we explore the kinetics and phase separation of ternary lipid mixtures of saturated lipid, unsaturated lipid, and cholesterol. Order parameters employed in the Cahn-Hilliard theory provide insight into the kinetics and mechanism of lipid phase separation. We observe three distinct time regimes in the phase separation process: a shorter time exponential phase followed by a power law phase followed by a longer time plateau phase. Comparison of lipid, protein and lipid-protein dynamics between all-atom and coarse-grained models identifies both quantitative and qualitative differences and similarities in the phase separation kinetics. Moreover, timescaling of dynamics of AA and CG simulation yields a similar kinetic mechanism of phase separation. The findings of this study elucidate fundamental aspects of membrane biophysics and the ongoing efforts to define the role of lipid rafts in the structure and function of cellular membrane.

Influenza A virus (IAV) infection relies on the action of the hemagglutinin (HA) and neuraminidase (NA) membrane proteins. The HA ligands anchor the IAV virion to the cells surface by binding the sialic acid (SA) present on the hosts receptors while NA is an enzyme capable of cleaving the SA from the extracellular environment. It is believed that the activity of NA ligands increases the motility of the virions favoring the propagation of the infection. In this work, we develop a numerical framework to study the dynamics of a virion moving across the cell surface for timescales much bigger than the typical ligand-receptor reaction times. We find that the rates controlling the ligand-receptor reactions and the maximal distance at which a pair of ligand-receptor molecules can interact greatly affect the motility of the virions. We also report on how different ways of organizing the two types of ligands on the virions surface result in different types of motion that we rationalize using general principles. In particular, we show how the emerging motility of the virion is less sensitive to the rate controlling the enzymatic activity when NA ligands are clustered. These results help to assess how variations in the biochemical properties of the ligand-receptor interactions (as observed across different IAV subtypes) affect the dynamics of the virions at the cell surface.

The genome contains several high-affinity non-functional binding sites for transcription factors (TFs) creating a hidden and unexplored layer of gene regulation. We investigate the role of such "decoy sites" in controlling noise (random fluctuations) in the level of a TF that is synthesized in stochastic bursts. Prior studies have assumed that decoy-bound TFs are protected from degradation, and in this case decoys function to buffer noise. Relaxing this assumption to consider arbitrary degradation rates for both bound/unbound TF states, we find rich noise behaviors. For low-affinity decoys, noise in the level of unbound TF always monotonically decreases to the Poisson limit with increasing decoy numbers. In contrast, for high affinity decoys, noise levels first increase with increasing decoy numbers, before decreasing back to the Poisson limit. Interestingly, while protection of bound TFs from degradation slows the time-scale of fluctuations in the unbound TF levels, decay of bounds TFs leads to faster fluctuations and smaller noise propagation to downstream target proteins. In summary, our analysis reveals stochastic dynamics emerging from nonspecific binding of TFs, and highlight the dual role of decoys as attenuators or amplifiers of gene expression noise depending on their binding affinity and stability of the bound TF.

We study the effects of osmotic shocks on lipid vesicles via coarse-grained molecular dynamics simulations by explicitly considering the solute in the system. We find that depending on their nature (hypo- or hypertonic) such shocks can lead to bursting events or engulfing of external material into inner compartments, among other morphology transformations. We characterize the dynamics of these processes and observe a separation of time scales between the osmotic shock absorption and the shape relaxation. Our work consequently provides an insight into the dynamics of compartmentalization in vesicular systems as a result of osmotic shocks, which can be of interest in the context of early proto-cell development and proto-cell compartmentalisation.

Binding reactions in effectively one-dimensional systems, such as proteins diffusing along DNA or other filaments, pose a fundamental coarse-graining challenge because stochastic trajectories are recurrent in one dimension and therefore do not admit a unique, separation-independent macroscopic association rate. As a result, continuum rate equations are not exact in 1D even for initially homogeneous systems. Here we develop a practical framework for mapping stochastic 1D reaction-diffusion dynamics onto effective kinetic models. Using mean-first-passage arguments and particle-based simulations, we define a density-dependent association rate and a corresponding single-rate approximation, and quantify when each provides an accurate description of the underlying stochastic dynamics. We implement 1D reaction-diffusion with excluded volume in the NERDSS software using a free-propagator reweighting algorithm and validate it against known pairwise and many-body limits. Our results show that ordinary rate equations with a single effective rate can accurately reproduce 1D reaction kinetics when the dimensionless parameter governing the ratio of intrinsic to diffusion-limited reactivity is small, with excellent agreement in the strongly rate-limited regime and increasing deviations as diffusion control strengthens. We further show that excluded volume in 1D can appreciably alter both kinetics and equilibrium populations, even at modest particle densities, by reducing accessible length and introducing blockade effects. Together, these results provide quantitative guidance for selecting between spatial simulations, density-dependent rate models, and single-rate continuum descriptions of reversible 1D binding reactions.

We present a Monte Carlo framework for simulating the nucleation of biomolecular condensates in the grand canonical ensemble, which overcomes the limitations of fixed particle number and allows direct control of the dilute-phase concentration. Our approach combines conformation sampling with bias-enhanced cluster size sampling, enabling accurate sampling of individual cluster sizes in various conformations. Our method resolves nucleation free energy surfaces and is capable of capturing both classical and nonclassical nucleation mechanisms. We validated our method by reproducing structural properties of disordered proteins across a diverse benchmark set and by applying it to study nucleation in two phase-separating proteins, FUS-LC and NDDX4, using both HPS and MPIPI coarse-grained force fields. While both proteins exhibit classical nucleation behavior under the HPS model, only FUS-LC remains classical with MPIPI. In contrast, NDDX4 shows a distinctly nonclassical nucleation pathway under MPIPI, characterized by a metastable intermediate state near the cluster size of 53 protein chains. Morphological analysis reveals that clusters up to this point are compact and spherical, whereas larger clusters adopt a nonspherical, two-lobed geometry indicative of frustrated coalescence. These findings underscore the critical role of sequence composition and force field parameterization in shaping nucleation pathways and demonstrate the utility of our framework for uncovering complex, mechanism-rich free energy landscapes in biomolecular condensation.

The effect of salt on coacervation of synthetic or biological polyelectrolytes and polyampholytes is well-studied. However, recent experiments showed that largely uncharged IDPs (like FUS) also undergo LLPS at physiological salt concentrations such as [Cion][~]0.15M, dissolve at higher salt concentration, and again phase separate at even higher salt concentrations such as, [Cion][~]3M. Here we use analytical theory and simulations to reveal the mechanism of these transitions. At low [Cion], the ionic solution acts as a highly correlated medium conferring long-range effective attractive interactions between spatially distant monomers. In this regime the ion concentration inside the condensate is higher than in the bulk solution. As [Cion] increases, the correlation length in the ionic plasma decreases, and the condensate dissolves. Second LLPS at high [Cion] is due to the entropy-driven crowding, and ion concentration inside the condensate is lower than in the bulk. Our study unravels a general physical mechanism of salt-dependent reentrant behavior in LLPS in uncharged IDPs.

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

1Structural changes in proteins allow them to exist in several conformations. Non-covalent interactions with ligands drive the structural changes, thereby allowing the protein to perform its biological function. Recent findings suggest that many proteins are always in an equilibrium of different conformations and that each of these conformations can be formed by both the ligand-free and ligand-bound protein. By using classical statistical mechanics, we derived the equilibrium probabilities of forming a conformation with and without ligand. We found, under certain conditions, that increasing the probability of forming a conformation by the ligand-free protein also increases the probability of forming the same conformation when the protein has a ligand bound. Further, we found that changes in the conformational equilibrium of the ligand-free protein can increase or decrease the affinity for the ligand.

Pulsed Dipolar ESR Spectroscopy (PDS) is a uniquely powerful technique to characterize the structural property of intrinsically disordered proteins (IDPs) and polymers and the conformational evolution of IDPs and polymers, e.g. during assembly, by offering the probability distribution of segment end-to-end distances. However, it is challenging to determine distance distribution P(r) of IDPs by PDS because of the uncertain and broad shape information that is intrinsic to the distance distribution of IDPs. We demonstrate here that the Srivastava-Freed Singular Value Decomposition (SF-SVD) point-wise mathematical inversion method along with wavelet denoising (WavPDS) can aid in obtaining reliable shapes for the distance distribution, P(r), for IDPs. We show that broad regions of P(r) as well as mixed narrow and broad features within the captured distance distribution range can be effectively resolved and differentiated without a priori knowledge. The advantage of SF-SVD and WavPDS is that the methods are transparent, requiring no adjustable parameters, the processing of the magnitude for the probability distribution is performed separately for each distance increment, and the outcome of the analysis is independent of the users judgement. We demonstrate the performance and present the application of WavPDS and SF-SVD on model ruler molecules, model polyethylene glycol polymers with end-to-end spin labeling, and IDPs with pairwise labeling spanning different segments of the protein tau to generate the transparent solutions to the P(r)s including their uncertainties and error analysis.